arXiv:1503.06585vl [physics.ao-ph] 23 Mar 2015 
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Abstract 

The controllability of advection-diffusion systems, subject to uncertain initial val¬ 
ues and emission rates, is estimated, given sparse and error affected observations of 
prognostic state variables. In predictive geophysical model systems, like atmospheric 
chemistry simulations, different parameter families influence the temporal evolution of 
the system. This renders initial-value-only optimisation by traditional data assimilation 
methods as insufficient. In this paper, a quantitative assessment method on validation 
of measurement configurations to optimize initial values and emission rates, and how 
to balance them, is introduced. In this theoretical approach, Kalman filter and smoother 
and their ensemble based versions are combined with a singular value decomposition, 
to evaluate the potential improvement associated with specific observational network 
configurations. Further, with the same singular vector analysis for the efficiency of 
observations, their sensitivity to model control can be identified by determining the 
direction and strength of maximum perturbation in a finite-time interval. 

Keywords: Atmospheric transport model, emission rate optimisation, observability, obser¬ 
vational network configuration, singular value decomposition, data assimilation 

1 Introduction 

Air quality and climate change are influenced by the fluxes of green house gases, reactive 
gas emissions and aerosols in the atmosphere. The ability to quantify variable, yet hardly 
observable emission rates is a key problem to be solved for the analysis of atmospheric 
systems, and typically addressed by elaborate and costly field campaigns or permanently 
operational observation networks. The temporal evolution of chemistry in the atmosphere 
is usually modelled by atmospheric chemistry transport models. Optimal simulations are 
based on techniques of combining numerical models with observations. In meteorological 
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forecast models, where initial values are insufficiently well known, while exerting a high 
influence on the model evolution, this procedure is termed data assimilation (J8l). There 
is no doubt that the optimization of the initial state is always of great importance for the 
improvement of predictive skill. However, especially for chemistry transport or greenhouse 
gas models with high dependence on the emissions in the troposphere, the optimization 
of initial state is no longer the only issue. The lack of ability to observe and estimate 
surface emission fluxes directly with necessary accuracy is a major roadblock, hampering 
the progress in predictive skills of climate and atmospheric chemistry models. In order to 
obtain the Best Linear Unbiased Estimation (BLUE) from the model with observations, ef¬ 
forts of optimization included the emission rates by spatio-temporal data assimilation have 
been made. The first full chemical implementation of the 4D-variational method for at¬ 
mospheric chemistry initial values is introduced in |[9]|. Further, Elbern et al. ( iflTft f took 
the strong constraint of the diurnal profile shape of emission rates such that their ampli¬ 
tudes and initial values are the only uncertainty to be optimized and then implemented it by 
4D-variational inversion. This strong constraint approach is reasonable because the diurnal 
evolution of emissions arc typically much better known than the absolute amount of daily 
emissions. Moreover, several data assimilation strategics were designed to adjust ozone ini¬ 
tial conditions and emission rates separately or jointly in (23l l. Bocquet et al. introduced a 
straightforward extension of the iterative ensemble Kalman smoother in Q. 

In many cases, the better estimations of both the initial state and emission rates are not 
always sustained based on appropriate observational network configurations when using 
popular data assimilation methods, such as 4D-variation and Kalman filter and smoother. It 
may hamper the optimization by unbalanced weights between the initial state and emission 
rates, which can, in practice, even result in degraded simulations beyond the time inter- 
vall with available observations. The ability to evaluate the suitability of an observational 
network to control chemical states and emission rates for its optimised designis the a key 
qualification, which needs to be adressed. 

Singular value decomposition (SVD) can help identifying the priorities of observations 
by detecting the fastest growing uncertainties. The targeted observations problem is an im¬ 
portant topic in the field of numerical weather prediction. Singular vector analysis based on 
SVD was firstly introduced to numerical weather prediction by Lorenz ( Il2ll ). who applied 
it to analyse the largest error growth rates in an idealised atmospheric model. Because of the 
high cost of computation, the singular vector analysis was not widely applied until 1980s. 
Later the method of singular vector analysis of states of the meteorological model with high 
dimension was feasible (@). 

In atmospheric chemistry, studies about the importance of observations are still sparse. 
Khattatov et al. ( 11710 firstly analysed the uncertainty of a chemical compositions. Liao et 
al. (ED) focused on the optimal placement of observation locations of the chemical trans¬ 
port model. However, singular vector analysis for atmospheric chemistry with emissions is 
different since emissions play an similarly important role in forecast accuracy with initial 
values. Goris and Elbern ( 1141 ) recently used the singular vector decomposition to deter¬ 
mine the sensitivity of the chemical composition to emissions and initial values for a variety 
of chemical scenarios and integration length. 

Hence, in this paper, applying the Kalman filter and smoother as the desirable data 
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assimilation method we introduce an approach to identify the sensitivities of a network 
to optimize emission rates and initial values independently and balanced prior to any data 
assimilation procedure. Through singular value decomposition and ensemble Kalman filter 
and smoother, the computational cost of this approach can be reduced so that it is feasible 
in practice. Then, by the equivalence between 4D variation and Kalman filter for linear 
models, the approach is also feasible for the data assimilation of adjoint models via 4D 
variational techniques. 

This paper is organized as follows. In section [2j we describe the atmospheric transport- 
diffusion model with emission rates first and then reconstruct the state vector such that 
the emission rates are included dynamically. In section [3j the theoretical approach derives 
in order to determine the efficiency of observations or observational network configurations 
before running any data assimilation procedure. In section|4j based on the theoretical analy¬ 
sis in section [3j we discuss the ensemble approach to evaluate the efficiency of observation 
configurations and present elementary examples. In section [5] we present the approach 
to identify the sensitivity of observations by determining the directions of maximum per¬ 
turbation growth to the initial perturbation. In the appendices, the above approaches are 
generalised to continuous-time systems for comprehensive applications. 

2 Model description 

The chemical tendency equation including emission rates, propagating forward in time, is 
usually described by the following atmospheric transport model 


dc ... . . 

it = A{c)+ e(t) 


where A is a nonlinear model operator, c(t) and e(t) are the state vector of chemical con¬ 
stituents and emission rates at time t, respectively . 

The a priori estimate of the state vector of concentrations c(t) is given and denoted by 
Cb(t), termed background state. The a priori estimate of emission rates are usually taken 
from emission inventories, denoted by e/,(f). 

Let A be the tangent linear operator of A, the evolution of the perturbation of states c(t) 
and e(t) follows the tangent linear model with A as 


ddc . r 
—— = A dc + 5e(t), 


( 1 ) 


where 5c(t) is the perturbation evolving from the perturbation of initial state of chemical 
state Sc(to) = c(t o) — q,(to) and emission rates Se(t) = e(t) — e/,(t). 

After discretizing the tangent linear model in space, let M(-, •) be the evolution operator 
or resolvent generated by A. It is straightforward to obtain the linear solution of ([B with 
continuous time as 



( 2 ) 
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where 5c(t) € JR", 5e(t) € IR", n the dimension of the partial phase space of concentrations 
and emission rates. Obviously, M(-, •) € IR nxn . 

In addition, let y(t) be the observation configuration of c(t) and define 

Sy(t) = y{t ) -H(t)c b (t), 

where T~L{t) is a nonlinear forward observation operator mapping the model space to the 
observation space. Then by linearising the nonlinear operator T~L as H, the linearised model 
equivalents of observation configurations can be presented as 


5y(t) = H(t)6c(t ) + u(t), 

where 5y(t) € R m W, m[t) the dimension of the phase space of observation configurations 
at time t. u(t) is the observation error at time t of the Gaussian distribution with zero mean 
and variance R(t) € R m ( i ) xm ( i ). 

It is feasible to apply the Kalman filter and smoother into the model without any ex¬ 
tension if the emission rates are accurate, which implies the initial state of concentration 
is the only parameter to be optimized. However, if the emission rates are poorly known, 
they should be combined into the state vector so that both of them can be updated by a 
smoother application. To establish the model with a new combination of the initial state and 
emissions, let us rewrite the background of emission rates into the dynamic form 


e h (t) = M e (t, s)e b (s), 


where e b (-) is a n-dimensional vector of which the i th element is denoted by e’ b (-) and 
M e (t, s) is the diagonal matrix defined as 


M e (t,s ) 


( e b (*) 

e U s ) 



V 


\ 


e £0) / 


Since emission rates follow the diurnal variation, by taking the diurnal profile of emis¬ 
sion rates as a constraint, the amplitude of emission rates can be estimated by constant 
emission factors f lfTTI t. We reconstruct the dynamic model of emission rate perturbation as 


5e(t) = M e (t,s)6e(s). 


Then, (O can be written as 


6c(t) 


M(t,t 0 )Sc(t 0 ) + 



M(t, s)M e (s , to)6e(to)ds. 


Hence, we obtain the extended model with emission rates 


( Sc(t) M(t,t 0 ) // 0 M(t, s)M e (s,to)ds \ ( 5c(t 0 ) \ 

V 5e ( f ) ) V 0 M e (t,t 0 ) ) V ) ' 


(3) 

(4) 
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Typically, there is no direct observation for emissions. Therefore, we reconstruct the 
observation mapping as 

5y(t) = [H(t),O nxn ] ^ ^ + K<)> 

where 0 nX n is a n x n matrix with zero elements. 

It is clear now that both concentrations and emission rates arc included into the state 
vector of the block extended model ([5]), such that the Kalman smoother in a fixed time 
interval [to, fjv] can be applied to optimize both of them. In our approach (0. the dynamic 
model of emission rates is forced to follow the background evolution of emission rates. 
In fact, it was stated in several studies (Q, l fT3l ) that the best linear unbiased estimation 
(BLUE) of a random variable x, which implies this estimation, can minimise its variance. 
The estimate via fix-interval Kalman smoother is the BLUE depending on all observation 
configurations in time interval [^o 5 ^tv]- In our case of emission rates, the estimation of e(t) 
by Kalman smoother on [to, Tv] can be represented generally as the conditional expectation 
E[e{t)\{y(t), t G [to - Tv] }]. By the linear property of conditional expectation, 

E[e(t)\{y(t),t € [t 0 ,t N }}} 

= E[M e (t,s)e(s)\{y(t),t G [fo,Uv]}] = M e (t,s)E[e(s)\{y(t),t € [fo,Uv]}], 

which implies the dynamic model of emission rates ([3]) satisfies the constraint of the diurnal 
shape of emission rates if [to, tjv] covers 24 hours. 

3 Efficiency of observation networks of atmospheric inverse mod¬ 
elling with emission rates 

As mentioned before, the observational network configurations cannot necessarily help im¬ 
proving the initial state and emission rates in a balanced way. If the estimation of both 
initial state and emission rates can be improved significantly, we call the corresponding 
observation configurations as efficient or of high efficiency for both. Otherwise, the ob¬ 
servation configurations are only efficient to initial state or emission rates. However, it is 
usually difficult to foresee the efficiency of observation configurations. Hence, the lack of 
the knowledge of the efficiency of observations may lead us to give the poor initial guesses 
and waste computational resource. In this section, we will introduce the theoretical ap¬ 
proach to determine the efficiency of observations via the Kalman filter and smoother in a 
finite-time interval. 

3.1 Theoretical analysis for the general discrete-time system 

For the application in atmospheric chemistry, let us consider the discrete-time system first. 

Generalizing the extended atmospheric transport model with emission rates in a discrete 
time internal [to, fi, ■ ■ ■ , tjv] to the following abstract linear system: 

x{t k+ 1 ) = M(t k+ i,t k )x(t k ) +e(t k ), (6) 
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y{t k ) = H(t k )x(t k ) + v(t k ) 


(V) 


where x(-) € R n is the state variable, y(t k ) € JT m( P 1 is the observation vector at time 
t k , the model error e(t k ) and the observation error v(t k ), k = 1, • • • , N follow Gaussian 
distribution with zero mean and Q(t k ) and R(t k ) are their covariance matrices respectively. 

Denote the estimation of x(t k ) based on {y(to), ■ ■ ■ , y(t k )} by x(t k \t k ), termed as the 
analysis estimation, the estimation of x(t k ) based on {y(to), ■ ■ ■ ,y(t k -i)} by x(t k \t k -i), 
termed as forecast estimation. Correspondingly, P(t k \t k ) and P(t k \t k -i) are the analysis 
error and forecast error covariance matrices of x(t k \t k ) and x(t k \t k -\) respectively. For 
convenience, the main results of the discrete-time Kalman filter can be summarised as fol¬ 
lows: 

(1) Analysis step: 

K(t k ) = P(4|4_i)iF T (4)(iF(4)P(4|4_i)iF T (4) + P(4)) -1 ; 
x{t k \t k ) = x(4|4_i) + K(t k )(y(t k ) - H(t k )x(t k \t k _i)); 

P(4|4) = (I — K{t k )H(t k ))P(tk\tk-i)\ 

(2) Forecasting step: 

x(t k+l \t k ) = M(t k+ i,t k )x(t k \t k ); 

P(t k+1 \t k ) = M(t k+ i,t k )P{t k \t k )M T (t k+ i,t k ) + Q{t k ), 

where for any matrix M, M T is the adjoint of M and AP 1 is the inverse of M. 

Denote the first guess of initial variance as P(fo|f-i) and select P(fo|f-i) and R(t k ) to 
be symmetric and positive definite. Then we can rewrite 

P(4|4) = P(4|4-i) 

- P(4|4-i)P T (4)(-ff(4)P(4|4-i)P T (4) + P(4)) _1 P(P)P(4|4-i), 

and by the matrix inverse lemma dl24l f. we have 

P _1 (4|4) = P -1 (4|4-i) + H(t k ) T R~ l (t k )H(t k ). (8) 

Further, assume the model error, which is usually unknown, is negligible. Then, we obtain 
P _1 (ffc+i|4) = M~ T {t k+ i,t k )P(t k \t k y 1 M~ l (t k+1 ,t k ). (9) 

ffence, by the deduction based on © and t©, we have 
P _1 (f fc+ i|4) 

= Af _T (4 + i, 4)P _1 (4|ffc-i)Af _1 (ffc + i,ffc) 

+M~ T (t k+ i,t k )H T (t k )R~ 1 (t k )H(t k )M~ 1 (t k+1 ,t k ) 

= M _T (4 + i,4_ 1 )P“ 1 (4_i|f fc _ 2 )AF _1 (4 + i,4_i) 

+M _T (4 + i,4_i)fT T (4_i)P _1 (4_ 1 )iT(4_i)M' 1 (4 + i,4_ 1 ) 

+M~ T (t k+1 ,t k )H T (t k )R~ 1 (t k )H(t k )M~ 1 (t k+ i,t k ) 
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= M T (t k+ i,t 0 )P l {t Q \t-i)M 1 (t k+ i,t 0 ) 

k 

+ ^ M - T (t k+1 ,t i )H T (t i )R- 1 (t i )H{U)M- 1 (t k+l ,t i ). 

i =0 

Define x(to\t k ) = E[x(to)\y(to ),..., y(t k )] and denote its covariance matrix as 

P(t 0 \t k ) = E[(x(t 0 ) - x(t 0 \t k ))(x(t 0 ) - x(t 0 \t k )) T ], 

which, according to the definition of x(to\t k ), is the covariance of the estimate of the state 
from the fixed-interval Kalman smoother. Then 


P _1 (*o|4) (10) 

= E[M~ 1 (t k+ i,t 0 )(x(t k+1 ) - x(t k +i\t k ))(x{to) - x(t k+1 \t k )) T M~ T (tk+iR.o)] -1 
= M T (4 + i,to)P _1 (tfe + i|4)M(tfc +1 ,to) 

k 

= P _1 (t 0 |t-i) + ^2M T (ti,t 0 )H T (ti)R~ 1 (ti)H(ti)M(t i ,to). 
i =0 

In particular, for k = N, taking the observations in the entire time interval into account, 
we have 


N 

P-\t 0 \t N ) = p-\ f 0 |f-i) + 5] MT (tu to^MR-HumtiWfato). (11) 

i =0 

It is clear that (flTT) includes all known information of the model with initial variance 
and observation configurations before any data assimilation procedure. At the same time, it 
is independent of any specific data and states. Actually, if we define 


Q = 


( H(t 0 )M(to,t 0 ) \ 
V H(tN)M(t N ,t 0 ) ) 


and 


/ R~\to) 


Rr 1 = 


R~Hh) 


R ) 


cm can be written as 


( 12 ) 


P-\t 0 \t N ) = P-^olM + Q T R~ X Q, (13) 

where Q T R~ 1 Q equals the observability Gramian with 1Z 1 from control theory, which is 
appropriate for how well states of a model can be inferred by the external observations. 
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Though d meets the demand to represent the covariance by all known information 
before starting the data assimilation procedure, the statistic interpretation of the inverse of 
covariance is still blurred to the application. Therefore, for evaluating the improvement of 
the estimation with the initial variance P(to|f-i), we define relative improvement covari¬ 
ance as 


P - P(t 0 \t N ))P 2(f 0 |f-l) 

= I - P~^{t 0 \t- 1 )P(t 0 \t N )P~^(t 0 \t-i), (14) 

where I is the identity matrix. 

The above improvement covariance is a normalised matrix of the difference between 
the initial variance P(to|t-i) and the covariance matrix P(fo|tiv) from Kalman smoother. 
Especially, P _ 2 (t 0 |t_ 1 )T , (fg|f 7 v)P ~2 (t 0 |i_ 1 ) can be understood as the covariance matrix 
from the fixed-interval Kalman smoother normalised by the initial variance. The symmetric 
normalised matrix guarantees the improvement covariance to be positive-definite. Further, 
its singular values or the eigenvalues are bounded since the sum of the eigenvalues of a 
matrix is equal to its trace. In fact, from (fl4l) . we have 

0 ^ P~^(t 0 \t-i)(P(t 0 \t-i) - P(t 0 \t N ))P~^(t 0 \t-i) < I , 

which implies that its trace is always less than the dimension of the state vector of the model. 

Since P(io|f rv) is unknown before the data assimilation procedure is finished, we rewrite 
the relative improvement covariance as 

P _ 5(i 0 |t_ 1 )(p(t 0 |t_ 1 ) - P(t 0 \t N ))P-^(t 0 \t_i) 

= p~^(t 0 \t-i)(p(t 0 \t-i) - (p _1 (t 0 |t- 1 ) + g t p~ 1 g)~ 1 )p~^{ t 0 \t-i) 

= i-p-5( io |t_ 1 )(P” 1 (f 0 |f-i) + e? T p- 1 a)' 1 P'5(i 0 | i _ 1 ) 

= i - + p^{to\t- 1 )G T n- l gp^{t 0 \t-i))- 1 . (is) 

It is clear to see from (fl5T ) the improvement covariance defined in (fl4l ) that 
I + P^ (iolt-i^ft-^P* (t 0 |t_i) 

is still invertible even without observability of the system, which means Q T Q is not full- 
rank. However, it is not easy to calculate (fl5T ) directly. Hence, by singular value decompo¬ 
sition, 

P^(fo|f-i)0 T P~^ = VSU T , 

where V and U are unitary matrices consisted of the left and right singular vectors, S is the 
rectangular diagonal matrix consisting of the singular values. 

We can simplify (fl5l) as 

P - ^(f 0 |f-i)(P(fo|t-i) - P~%(t 0 \t N ))P(to\t_i) 

= i — (i + p^ (t 0 \t^i)g T iz- 1 gp^ (toit-i))- 1 
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= I — (/ + VSS T V T )~ 1 

= vv T -ivv T + vss T v T )~ 1 
= vv T -ivii + ss 7 ^)- 1 

= V{I~{I + SS T )~ 1 )V T 


- E 


1 + s 


2 v i v i 


(16) 


where r is the rank of (fl5l) and v h is the i th left singular vector in V related to the singular 
value Si, which is the i th element on the diagonal of S. 

Let us consider the improvement of each element in the state vector as the corresponding 
value in the diagonal of the relative improvement covariance. From (fl6l) . we denote the 
relative improvement of j th element in a:(to) as Pj, then, 


S-E 



where Vij is the j th element of v % . 

Get a deeper insight into the capacity of the observation networks to improve the esti¬ 
mation of all states of the model, some important indices need to be considered. In fact, 
in order to evaluate the total improvement of the model, the nuclear norm for matrices, or 
equivalently, the 1-norm is appropriate, which is defined as 

||M||i = tr (VM T M), 


where M is any matrix and tr(-) denote the trace of the matrix. 

For (fT6l) . we denote 

P = p-h(t 0 \t-i)(P(t 0 \t-i) - P(t 0 \t N ))P~^ (t 0 \t-i), (17) 


according to (fl6l) . 


|P||i = 


E 

i =1 


1 + S. 


2 > 


which is called the total improvement value. 

As we mentioned before, ||P||i < \\I\\i = n, where n can be considered as the total 
improvement value, if the system is fully known, which implies the optimal estimation is 
the value of state. Thus, if we consider the ratio 


P = 




€ [ 0 , 1 ), 


(18) 


the percentage of the total improvement of the model is obtained, which is called the relative 
improvement degree. 
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3.2 Application to the extended atmospheric transport model with emissions 

For the atmospheric transport model extended with emissions, composing the dimension of 
the original state c E R” and emission rates e E R n respectively, we divide (fT71) into a 
block matrix 


P = 



pee 

pe 


2 n 

E 

2=1 


1 + sj 


{v\ 


€ R 


2nx2n 


where P c is the relative improvement covariance of the state c(to)> P e is the relative im¬ 
provement covariance of the emission rates e(fo), P ce = (P ec ) 1 is the relative improvement 
covariance between c(fo) and e(to) and (vf , vf ) T = v^. 

Then, it is easy to calculate 


P c = 


2 n 

E 

2=1 


1 + 5 


—vfivP 
2 2 2 


P C = 


2 n 

Sr 


2=1 


+ 5 


2 u l 2 


Further, the relative improvements of j th element in c(to) and e(fo) are 


P i = 


2 n 


Sr 


2=1 


+ sf 


to 




p e = 

3 


2 n 


+ sf 


tofc) 


where and v\ rj are respectively the j th element of vf and vf. 

Moreover, the total improvement values of concentration and emission rates are 


|R c ||i = 


2 n 


. , 1 + ■ 
2=1 




l^lli = 


2 n 

Sr 

2=1 


+ 5; 


rtrtofwf 


It is worth noting that 

P c = (P c (t 0 \t^))-HP c (t 0 \t-i) - P c (2 0 |tiv))( J P c (to|t-i))-^ 


and 

P e = (P e (io|t-i))-^(R e (to|t-i) - P e {to\tN)){P e (t 0 \t- 1 ))-$ 

if and only if there is no correlation between the initial concentration and emission rates. In 
fact, if we assume P ce (fo|^-i) = 0 nxn , the coiTesponding relative improvement degrees of 
concentration and emission rates are defined as 



n 


P 



n 


According to (flSl) . it is obvious that ff E [0,1) and [/' E [0,1) show the percentages 
of the relative improvements of concentration and emission rates, respectively. However, 
since 



n n 


> 1, 
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which indicates the normalisation is just with respect to the extended covariance matrix 
rather than specified to the state c and emission rates e. The relative improvement degree 
cannot serve our objective to distinguish the observability of concentration and emission 
rates and balance them quantitatively. However, by observing the block form of P, it is 
easy to obtain 

ll^li +ll-P e ||i = ll^lli- 

For comparing the improvement of the concentration and emission rates, we define relative 
improvement ratios for the state or the emission rates as 


I pc 


I P e 


p = 


\p\ 


p = 


\p\ 


p e + p c = 1. 


In a sum, if the total improvement value or relative improvement degree of the model is 
almost zero, the relative improvement ratios do not need to be considered since no state of 
the model is improved. Otherwise, {P?}” =1 and {PJ}” =1 , which show the improvement of 
each parameter j of concentrations and emission rates respectively, can help us determining 
which parameters can be optimized by the existing observation configurations. By compar¬ 
ing p c with p e , it is clear that the estimation of the one with the larger ratio (larger magnitude 
of the improvement covariance) can be improved more efficiently by the existing observa¬ 
tion configurations. In other words, if p c > p e , the existing observation configurations are 
more sensible to the initial values of concentration. Conversely, if ff < p e , the observation 
configurations can help improving the estimation of emission rates more. According to p c 
and p e , the ’weight’ between the concentrations and emission rates can be decided quanti¬ 
tatively. In a data assimilation context, where observations are in a weighted relation to the 
background, the BLUE favours the more sensitive parameters. 

In a special case that p e is very close to zero, the emission rates can be viewed as the 
input in the model without any optimization when the data assimilation procedure is started. 

For the completion of the theorem and wider application, the generalisation of the above 
method for the continuous-time system is introduced in Appendix A. While the derivation is 
totally different from the discrete-time system, it can be shown that the theoretical analysis 
still holds for the continuous-time system. 


4 Application to the ensemble Kalman filter and smoother 

In practice, the standard Kalman filter and smoother cannot be applied directly to transport 
modells due to their computational complexity. The ensemble Kalman filter (EnKF) and 
smoother (EnKS), as a Monte Carlo implementation originating from Kalman filter and 
smoother, are suitable for problems with a large number of control variables and are an 
important tool in the field of data assimilation ( ifTHl f. In this section, we will introduce how 
to apply the above method if an Ensemble Kalman filer and smoother are applied as the data 
assimilation method. 

For the abstract discrete-time system ([6]), we denote the ensemble samples of |L-i ) 
and x{ti\ti) , i = 1, ■ • • , N respectively by 

X(ti\ti-l) = l),X 2 (L|fi—l), • • • ,X q {ti\ti-l)), 
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X(ti\ti) = (xi(ti\ti),X2(ti\ti),--- ,X q (ti\ti)), 

where q is the number of ensemble members. 

Correspondingly, their ensemble means arc 

1 q 1 

x{ti\ti— i) — — ^ ' X k (ti\ti— i) — — X (ti |f i—i )lqx 1 5 

1 q 1 

x(ti\ti) — — ^ ' X k (ti |f j) — —A (f j|fj)lgxl5 

q tl q 

where t lX] is a i x j matrix where each element is equal to 1 . 

Note the ensemble perturbation matrix consist of the perturbation of each sampling by 

X(ti\U-i) = X(ti\ti-i) - -X(ti\ti-i)\ qxq . 

q 

Then the ensemble covariance is 

= -^^X{ti\ti-i)X T P{ti\ti) = -^-jX(ti\ti)X T (ti\ti). (19) 

Further, we define the ensemble observation equivalent in the entire time interval in 
observation space as 

£ 

V k = Gx k (t 0 \t-i), k= 1, • • • ,q 

and denote the ensemble mean and the forecast error covariance matrix in observation space 
by 

y f = 152 y£> P L = ~i - y f ^yl ~ y f "> T = sp(to\t-i)g T . 

q k=1 q k=1 

Similarly, we denote the ensemble covariance between the initial states and the forecasting 
observations by 

1 q 

p £ y = —y - x(t 0 \t-i))(yl - y f ) T = P(f 0 |f-i)£ r 

q 1 k =1 

In addition, we define the ensemble observations as 

Vk{U) = y(U ) + u k {U), k = 1, • • • , q, i = !,■■■ ,N 

and assume v(ti ) = | Y£k=i ,y fc(C) = 0. Denote the ensemble covariance of observation 
errors by R(ti) = ^ Yxk =l v k{U)v k {U)- Further, we assume 

R-\t 0 ) 

R-'ih) 

R-\t N ) 
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It is shown by Evensen ( lfl2l ') that the ensemble forecasting and analysis covariances 
have the same form with the covariances in the standard Kalman filter. It indicates that 
© and © are also true for P(ti\ti) and P{ti\ti-i). So upon substituting P(fo|f-i) into 
P(£o|t-i) in ( f1~3l ). according to matrix inversion lemma, we obtain 

P(to\tN) 

= p(f 0 |f_t) - + K-hP(to\t^ 1 )g T K-^)- 1 'ji- 1 igP(to\t- 1 ) 

= p(fo|t-t) - P£ y n~Hi + n- 1 2pf y n- 1 2)- 1 n- 1 Hpf y ) T . (20) 


Then, 

P~^(t 0 \t-i)(P(t 0 \t-i) - P(t 0 \t N ))P~^{t 0 \t- 1 ) 

= + (21) 

To simplify (I2TI) . let YliLi m (f«) = m, by singular value decomposition, we have 

p-z(t 0 \t- 1 )pf y K~* = VSU T € R nxm (22) 

where U € K, mxm consists of the eigenvectors of 'Rr^QP{to\t_\)Q T 'Rr^, V € R nxri 
consists of the eigenvectors of P2(t 0 |t_ 1 )^ / 7?. _1 ^P2(fo|f-i), S € IT nxm consists of the 
singular values. 

Let r be the rank of (1221) . the ensemble relative improvement covariance is defined as 


P 2 (to\t-i){P(t 0 \t-i) — P(to\tN))P 2 (to\t-[ 

VS T U T (UU T + [/(S5 t )[/ t )- 1 (/SI/ t 


T c \-l 


VS 1 (/ + S 1 S) 

r 2 

Si 


sv - 1 


i= 1 


+ 5 . 


2 v i v I 


(23) 


and its diagonal shows ensemble relative improvements of states. 

Obviously, (1231) has the same form as (fl6l) . In fact, for the linear dynamic model, we 
have 

p-kto\t-i)P£ y n- 12 = P^ (to\t-i)Q T P~^, (24) 

which implies if we just substitute 1Z and P(to\t-i) in (1T61) by P and P(io|t-i), the final 
results of (fl6l) and (1231) arc equivalent. However, it is much more efficient to calculate 
practically the singular values and vectors of the left-side matrix of (f24l) than to calculate 
the singular values and vectors of the right-side matrix of (l24l) . since the explicit calculation 
of Q is not necessary. 

Further on, it is worth noticing that if the ensemble size q is less than the dimension 
of the model n, the initial ensemble covariance P(f 0 |f-i) is not invertible. In this case, 
it is reasonable to replace P ~2 (to\t-i) by the pseudo inverse of P^(to\t-\), denoted by 
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pt 2 (t 0 |t-i), and calculate it by singular value decomposition. In fact, according to (fl9l) . 
we have 

P(fo|t-i) = ^-X(fo|t-i)X T (fo|t-i). (25) 

q - i 

Then, by singular value decomposition, 

—Uj=A-(f 0 |t-i) = VqSqUq , (26) 

where Vo E IR nxn and (To E lii qxq consist of the left and right singular vectors, So E R nx ' 1 ' 
is a rectangular diagonal matrix with singular values {sojl^oi ^ 0}f =1 on the diagonal. 
Thus, 

P(to\t-!) = VqSqUq UqSq V 0 t = VoSoSfiV 0 T = V 0 S 2 0 V 0 T , (27) 

where Sq = SqSq G R nxn is a block diagonal matrix with the rank ro and the diagonal 

(s 0 ir** ’ ^Oro’ x ( n—r o))* Hence, 

P^(t 0 \t- 1 ) = V 0 SlvT, 

where Sq is the pseudo inverse of So with the diagonal (l/soi,- -- , l/sor 0 >Oix (n-r 0 ))- 
Then, the ensemble relative improvement covariance can be rewritten as 

pH(f 0 |t-i)(P(to|i-i) - 

= P^(to\t-i)P(t 0 \t-i)P^(t 0 \t-i) - P^(to\t-i)P(to\t N )P^(t 0 \t-i) 

= Vb^JVB r (VbSgVo r )Vb^JvS r — J P t 5(to|t-i)^(to|tjv)^P + 3(to|t-i) 

= V 0 I ro V 0 T - P^(to\t-i)P(to\t N )P^-(to\t-i), 


where I ro is the diagonal matrix with the diagonal (lixr 0 > Oix(n-r 0 ))- 

It is clear from (l20l) that P f 2 (i 0 |t_ 1 )P(t 0 |^ 7V )P t 2 is still nonnegative definite 

while P(to\t-i) is not with full rank, so if we use the same notation of standard Kalman 
filter and smoother to denote the ensemble relative improvement covariance, which means 

P = P%o|f_i)(P(fo|t-t) - P(t 0 \t N ))P^(to\t-i), 


then, 0 n xn P: P < I ro . Further, the ensemble relative improvement degree is 


IPII 


IPI 


P = 


IT 


ro 111 


?"0 


- € [ 0 , 1 ). 


(28) 


As to the atmospheric transport model extended with emissions, for the distinction of 
the improvements for concentrations and emission rates, the ensemble relative ratios are 
still 

_ c || P c || i _ e ||P e ||i 

P = , P = 


P 


P 
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If we further consider the nonlinear dynamic model, we can renew the definition of the 
forecasting observation configurations as 


y{ = Q{x[(t 0 )), k =!,-■■ ,q, 


such that it can fully follow the nonlinear model, where Q is a nonlinear operator. 
Correspondingly, its ensemble mean and covariance are 



Thus, for a nonlinear dynamic model, the extended ensemble Kalman filter to the theo¬ 
retical approach in Section [3j the only approximation is the limited size of the ensemble. 

Hence, for the extended atmospheric transport model with emission rates, the analysis 
is similarto the analysis in Section 13.21 

4.1 Example 

Consider a linear advection-diffusion model with periodic horizontal boundary condition 
and Neumann boundary condition in the vertical direction on the domain [0,14] x [0,14] x 



where Sc, Se and Sd are the perturbations of the concentration, the emission rate and de¬ 
position rate of a species respectively. v x and v y are constants and K(z) is a differentiable 
function of height z. 

Assume At = 0.5, the numerical solution is based on the symmetric operator splitting 
technique ( |[25l0 with the following operator sequence 


5c(t + At) = T x TyD z AD z T y T x 5c{t) 


where T x and T y are transport operators in horizontal directions (x, y ), D z is the diffusion 
operator in vertical direction ( 2 ). The parameters of emission and deposition rates are 
included in A. The Lax-Wendroff algorithm is chosen as the discretization method for 
horizontal advection with Ax = Ay = 1. The vertical diffusion is discretized by Crank- 
Nicolson discretisation with the Thomas algorithm as solver. The horizontal domain is 
[0,14] x [0,14] with the horizontal space discretization interval, while the vertical domain 
is [0,4] with A z = 1. So the number of the grid points N g = N x X N y x N z = 1125, 
where N x = 15, N y = 15, N z = 5. 

In addition, we choose M e (t,to), a continuous function in time t, to formulate the 
temporal background evolution profile shape of the emission rate as 


e b {t) = M e (t,t 0 )e b (t 0 ) 


where e b (to) is the initial value of emission rate. 
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With the same assumptions of At and grid points in the 3D domain, the discrete dy¬ 
namic model of emission rates is 

5e(t + At, i,j, l) = M e (t + At, t)Se(t, i,j, l ), 
where l),i,j G {0, • • • , 14}, l G {0, • • • , 4}} arc the coordinates of grid points and 

M e (t + A t,t) = eb{t + A t)/eb{t). 

For expository reasons the background assumption of Sd is denoted by 5db, which is 
kept fixed. 

According to the discretization of the phase space, we always assume there is only one 
fixed observation configuration in this example. It indicates that the observation operator 
mapping the state space to the observation space is a 1 x 2 N g time-invariant matrix. 

Set 500 (the ensemble number N) samplings for the initial concentration and emission 
rate respectively by pseudo independent random numbers and make the states correlated by 
moving average technique. 

Advection test: For the advection test (Fig. 1 to Fig. 6), we assume the model with a 
weak diffusion process ( K(z ) = 0.5e _z ) and there is one single observation configuration 
of the concentration in the lowest layer at each time step, denoted by ’Obs-cfg of cone’ in 
figures. Besides, the emission source is assumed mainly from the location shown by the 
blue point in figures, named ’Emss-source’. 

If we set the data assimilation window to lOAf and the wind is from southwest, the 
left-side subplot in Fig. 1 shows the estimation of the concentration is probably improved at 
the field around the observation under the small assimilation window. Meanwhile, though 
the right-side subplot in Fig. 1 shows hardly improvement of the emission rate, we can see 
from the first line of Table 1 is feasible only for the concentration, for the simple reason that 
the single observation configuration cannot detect the emission within the corresponding 
assimilation window. 

If we consider the same case as Fig. 1, but now extending the data assimilation window 
to 35 A t. Fig. 2 shows the field where the concentration is potentially improved is enlarged 
since the states arc more correlated with the extension of assimilation window and the es¬ 
timation of the emission surrounding the emission source is improved, compared to the 
Fig. 1. The quantitative balance between the concentration and the emission is shown by 
the relative improvement ratios in the second line of Table 1. 

If we further extend the data assimilation window to 48Af, it is clear to see from Fig. 3 
and the third line of Table 1 that the states arc more correlated such that more areas can be 
analysed and improved by the single observation configuration. Meanwhile, the improve¬ 
ment of the emission is dominant with increasing time. 

Fig. 4 to Fig. 6 show the relative improvements of the concentration and emission rate, 
when the model domain is under a northeasterly wind regime, and assimilation windows 
of with the data assimilation windows lOAf, 35 At and 48Af respectively. It is easy to 
imagine that with northeasterly winds, whatever the duration of the assimilation window is, 
the emission is not detectable and improveable by the single observation configuration. This 


16 





Figure 1: Advection test with data assimilation (DA) window 10 At and southwest wind. 
Isopleths of ensemble relative improvements of the concentration and emission rate are 
shown in the leftside and rightside figures respectively. 


hypothesis is successfully tested by our approach, the results of which are clearly visible in 
Fig. 4 to Fig. 6 and Table 2. 

Emission signal test: The purpose of emission signal test (Fig. 7 and Fig. 8) is to show 
the approach is also sensitive to the different background profile of the emission rate evo¬ 
lution. Hence, the only distinction between the situations in Fig. 7 and Fig. 8 is the back¬ 
ground profile of the emission rate during the assimilation window 48A t. Actually, Fig. 8 
is the same case as Fig. 3. Thus, the result of the approach is clearly shown in Table 3 
that the strong emission signal or the distinct variation of the emission rate during the data 
assimilation window is significant to the model to recognize the source of the changes of 
the concentration and improve the estimation of the states. 

Diffusion test: The diffusion test (Fig. 9 and Fig. 10) aims to test the approach via 
comparing the ensemble relative improvements of the concentration and the emission rate 
of the model with a weak diffusion process and a strong diffusion process. For the case in 
Fig. 9, all assumptions are same with the situation in Fig. 2 except that the single observation 
configuration is at the top layer instead. The only difference of the assumptions between 
Fig. 9 and Fig. 10 is that K(z) = 0.5e -;z2 in Fig. 9 and K(z) = 0.5e -22 + 1 in Fig. 10. 

Comparing Fig. 2 with Fig. 9, it is obvious that the different observation location in¬ 
fluence on the distribution of the relative improvements of the concentration greatly. From 
Table 5, the total improvement value of the concentration in the lowest layer for Fig. 2 is 
shown to be larger than the one for Fig. 9. Besides, it can be seen in Table 4 that the ob¬ 
servation configuration in the top layer cannot detect the emission with such weak diffusion 
under the assimilation window 35A. 

If we compare Fig. 9 with Fig. 10, it is shown in Table 5 that both the total improvement 
value of the concentration in the lowest layer for Fig. 10 and the weight of the emission 
rate increase, which implies that the observation configuration is more efficient to detect the 
emission and improve the estimation of the state of the model with the strong diffusion in 
Fig. 10. 
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Figure 2: Advection test with DA window 35At and southwest wind. Plotting conventions 
are as in Fig. 1. 
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Figure 3: Advection test with DA window 48At and southwest wind. Plotting conventions 
arc as in Fig. 1. 




Figure 4: Advection test with DA window 10 At and northeast wind. Plotting conventions 
are as in Fig. 1. 
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Figure 5: Advection test with DA window 35 At and northeast wind. Plotting conventions 
are as in Fig. 1. 



Figure 6: Advection test with DA window 
are as in Fig. 1. 



x with v = -0.5 

X 

and northeast wind. Plotting conventions 
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Emss profile shape 





Figure 7: Emission signal test (weak) with DA window 35 At and southwest wind. Plotting 
conventions are as in Fig. 1. 


Emss profile shape 





Figure 8: Emission signal test (strong) with DA window 35 At and southwest wind. Plotting 
conventions are as in Fig. 1. 



20 












































Conc-imp under 35 time steps Emss-imp under 35 time steps 



Figure 9: Diffusion test (weak) with DA window 35At and southwest wind. Plotting con¬ 
ventions are as in Fig. 1. 


Conc-imp under 35 time steps Emss-imp under 35 time steps 



Figure 10: Diffusion test (strong) with DA window 35At and southwest wind. Plotting 
conventions are as in Fig. 1. 
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P c 

P e 

Fig. 1 

0.9979 

0.0021 

Fig. 2 

0.5107 

0.4893 

Fig. 3 

0.0345 

0.9655 


Table 1: 



P c 

P e 

Fig. 7 

0.6811 

0.3189 

Fig. 8 

0.5107 

0.4893 



P c 

P e 

Fig. 4 

0.9977 

0.0023 

Fig. 5 

0.9974 

0.0026 

Fig. 6 

0.9974 

0.0026 


Table 2: 



P c 

P e 

Fig. 9 

0.9977 

0.0023 

Fig. 10 

0.7755 

0.2245 


Table 3: 


Table 4: 


5 Sensitivity of observation networks to the initial state and emis¬ 
sion rates 

From the above discussion, we can determine the efficiency of the observation network 
by evaluating the improvement of estimation of initial state and emission rates separately, 
before we run the data assimilation by Kalman filer and smoother. However, it does not 
provide the information about the improved configurations of observations which can help 
improving the estimations. In this section, independent of any concrete data assimilation 
method, we will introduce the singular vector approach to identify the sensitive directions 
of observations to the initial state and emission rates. 

Consider the generalized discrete-time linear system: 


5x(t k+1 ) = M(t k+1 ,t k )5x(t k ), 

where 5x(to) = x(to) — x(to), x(to) is any estimate of x(fo)- 

Assume the observation mapping is accurate, which implies the data is the only source 
of observation errors, we have 


M4) = H(t k )5x{t k ). 

Define the magnitude of the perturbation of the initial state by the norm in the state 
space with respect to a positive definite matrix Wq 

||<5x(t 0 )||wo = (4'x(fo),W / ofe(f 0 ))- 

Similarly, we define the magnitude of the related observations perturbation in the time 
interval [to, • • • , tjv] by the norm with respect to a sequence of positive definite matrix 

{W(t k )}” =1 

N 

W s y\\{w(t k )} = ^2(^y(tk),w{t k )dy(t k )), 

k =0 

where Sy(t k ) = H(t k )5x{t k ). 
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pc 

low 

P e 

1 low 

Fig. 2 

2.8435 x 10" e 

2.3530 x 10' e 

Fig. 9 

2.3850 x 10 -7 

2.1627 x 10 -8 

Fig. 10 

7.8820 x 10 -7 

1.5946 x 10' 6 


Table 5: P[ 0V! and P[ olr are respectively the total improvement values of the concentration 
and emission rates in the lowest layer 


In order to find the direction of observation configuration which can minimize the per¬ 
turbation of the initial states, the ratio 


\\Sx(t 0 )\\ 2 Wo 

IMUw' 


5y ^0 


should be minimized. It is equivalent to maximize the ratio of the magnitude of observation 
perturbation and the initial perturbation 


\\ Sx (to)\\wo ’ 


Sx(t 0 ) / 0. 


Thus, we define the measure the perturbation growth as 


2 = H^ll \w{t k )} 

!l P*(*o)|lwb 

y ( Sy(t k ),W(t k )5y(t k )) 

{5x(t 0 ),W 0 5x(t 0 )) 

y (H(t k )5x(t k ),W(t k )H(t k )5x(t k )) 

(fe(fo),Wo<5x(to)) 

y (■ 6x(t k ),H(t k ) T W(t k )H(t k )6x(t k )) 
g {Sx(t 0 ), Wo<5x(f 0 )) 

y (5x{t 0 ),M(t k ,t 0 ) J H(t k )W(t k )H(t. k )M(t k ,t 0 )8x(to)) 

g (Sx(to),W 0 Sx(t 0 )) 

(5x(t 0 ),J2 k=0 M(t k ,t 0 ) T H(t k ) T W(t k )H(t k )M(t k ,t 0 )5x(t 0 )) 
{Sx(to),W 0 8x(t 0 )) 


(29) 


According to Liao and Sandu ( lfT8l l. singular vectors refer to the directions of the error 
growth in a descend sequence with respect to the descent singular values. Hence, in order 
to search the maximal directions of 


g 


2 


(<5x(t 0 ), G T WG5x(t 0 )) 
(5x(t 0 ),W 0 5x(t 0 )) 


8x(t 0 ) / 0, 
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where W = diag (W (to), ■ ■ ■ ,W ( tw )), Q and 7Z 1 have the same definitions with those in 
Section [3j we need to find out the solutions of the singular value problem: 

wpg T wgwp Vk = s\v k , gw 0 G T u k = s 2 k u k , 


where si ^ S 2 ^ ^ s n ^ 0, {v k }™ =l and {u k }f =1 are the corresponding orthogonal 

singular vectors. Then, max^( to w 0 9 2 = s '(- 

Especially, if the perturbation norms are provided by the choice Wq = P" 1 (fo|^-i) an d 
W = 1Z~ 1 defined in Section[3j 
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2 


(5x(to),g J n l Q5x(to)) 

(Sx(t 0 ),P- l (t 0 \t- 1 )Sx(to)y [0) * 


We need to search the directions of 


Pz(to\t-i)G T K 1 gp 1 2(to\t-i)v k = S 2 k v k ] (30) 

gP~ l (to\t-i)g T u k = s k u k , k = !,■■■ , n. 

Associated with (fl6l) . it is easy to find that the singular vector v k in (l30l) . which is 
the direction of k th -fast growth of the perturbation of observations evolved from the initial 
perturbation, is also the k th direction which maximize the improvement of estimation by 

Kalman filter and smoother (fl6l) . though the exact value of the eigenvalue of (fl6l) related to 

2 

v k is , rather than the eigenvalue s\ of (l30l) . Meanwhile, we can find that the leading 
singular value si is related to the operator norm of P as 

s 2 

IIPII = max llPxll = -~ o. 

||*||=1 1 + sf 

In addition, the similar analysis for the continuous-time system is presented in the appendix 

B. 


6 Discussion 

In the present work, approaches for determining the efficiency and sensitivity of observation 
configurations for the initial state and emission rates are established. Actually, to deal with 
the specific questions in atmospheric chemistry, some special operators are usually applied. 
For example, in order to consider the efficiency and sensitivity of observations in some 
certain locations, the local projection operator introduced by Buizza et al. (0|) can be 
applied into approaches in Section|4]and Section[5j 
Let L be the 0 — 1 diagonal matrix defined as 

I _ r 1 ) li € Lai 

11 ^0, otherwise. 

where L a is a fixed area and l t is the coordinate of i lh grid point. 
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To test the efficiency and sensitivity of observation configurations in a special area, by 
rearranging the observations y according to the locations, Q in (fT2l) should be defined as 


G 


( LH(to)M(to,to) \ 
LH(ti)M(ti,to) 


\ LH(t]y)M(tN,to ) / 


(31) 


If LH(-) is considered as the observation mapping, approaches in Section [3] and [5] can be 
applied. 

In addition, if there is a multiplication of emission rates in the following model 


-jj- = A Sc + B(t)Se(t ), 

then all approaches can also be applied into the extended model 

f 5c{t) \ = f M(t,t 0 ) f* o M(t,s)B(s)M e (s,t 0 )ds 

V Ht) ) V 0 Me (Mo) 


5c(t 0 ) \ 
5e(t 0 ) ) ' 


A The efficiency of observation networks for continuous-time 

systems 

Consider the abstract continuous-time system 

x{t) = M(t, t 0 )x(t 0 ) + s(t), 
y(t) = H(t)x(t ) + v(t), 

where x E !R, n is the state variable, y E R” 1 is observation vector at time t, the model 
error e(t) and the observation eiTor u(t), t E [to, M] follow Gaussian distribution with zero 
mean, while Q(t) and H(t) are their covariance matrices respectively. 

As in Section [3} we ignore the model error. It is well known that for the continuous 
Kalman filter, the covariance P(t) of the optimal estimation of the state at time t satisfies 
the integral Riccati equation 

P{t\t) = M K (t,to)P(to\t-i)Mx(t,to) + [ M K (t,s)K(s)R(s)K T (s)M'f < (t,s)ds, 

■n 0 

where K(t) = P(t\t)H(t)RT l (t) and Mx(t, to) = M(t, to)—fl Q M(t, s)K(s)H(s)Mk{s , to)ds. 

On one hand, 

M K (t,to)P{to\t-i)Mf < (t,to) 

= M(t,to)P(t 0 \t-i)MK(t,t 0 ) - [ M(t,s)K(s)H(s)M K (s,t 0 )P(to\t-i)Mf < (t,to)ds 

Jt 0 

= M(t,t 0 )P{to\t-i)Mj)(t,to) - f M(t,s)K(s)H(s)P(s\s)M^(t,s)ds 

Jt 0 
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Jto Jto 


On the other hand, 


'to 


M K (t, s)K{s)R{s)K T (s)M^(t, s)ds 


'to 


[. M(t,s ) — / M(t,rj)K(ri)H(r])MK('r],s)dri\K(s)R(s)K T (s)M^(t,s)ds 


M(t, s)K(s)R{s)K t (s)M%(t, s)ds 


’to 


M(t, ri)K(ri)H( t])Mk (rji s)K(s)R(s)K T (s)M^(t, s)dsdrj 
[ M (t, s)K(s)R{s)K T (s)M%(t, s)ds 

to 

M(t , s)K(s)H(s)Mk(s , rj)K(r])R(r])K T rj)dr]ds. 

Therefore, P(t\t) = M(t,to)P(to\t-i)M^(t,to). 

Since 




M 1 (t,t 0 ) = M K 1 (t,t 0 )M K (t,t 0 )M 1 (t,t 0 ) 

= M ^ 1 (t,to)[M(t,to) — [ M K (t,s)K(s)H(s)M(s,t 0 )ds\M- 1 (t,t 0 ) 

Jt 0 

= M^(Mo)- [ M^(s,t 0 )L(s)H(s)M{t,s)ds , 

Jt 0 

we obtain M^(t, to) = M _1 (t, to) + J^ q M ^ 1 (s,to)K(s)H(s)M ~ 1 (t, s)ds. 
Define x(to|t) = -E[x(to)|y(s), s € [to, t]] and denote its covariance matrix as 

P(f 0 |t) = £[(x(t 0 ) - z(t 0 |i))(x(t 0 ) - x(to\t)) T ], 


Hence, 


P~\t 0 \t) 

= [M~ l (t, t 0 )P(t\t)M~ T (t, t 0 )] _1 
= [P(t 0 |t-i)M^(t, t 0 )M~ T (t, t 0 )] -1 

= M T (t,t 0 )[M~ T (t,t 0 ) + [ M~ T (t,s)H T (s)K T (s)M^ T (s,t 0 )ds]P(t 0 \t-iy 1 

Jt 0 

= H _1 (to|t-i)+ [ M T (s,t 0 )H T (s)R~ 1 (s)H(s)M(s,t 0 )ds. 
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Let t = tw and define the observability mapping Q : R" —» £ 2 ([to, Lv]; R T ") as 

gf:=H(-)M(;t 0 )f, /er, 


its adjoint operator Q* is 

G*f = ~ f° M T (s,t 0 )H T (s)f(s)ds, f € L 2 ([t 0 ,t N ]-R m ). 

Jtjsi 

Further, we define PR 1 : £ 2 ([to, Lv]; R m ) —> £ 2 ([fi), Lv]; K™). 

^“7 := f € £ 2 ([t 0 ,tjv];K. m ). 

Thus, 

P~\t 0 \t N ) = + £*^“7, (32) 

where Q*TZ~ 1 Q is the observability Gramian of continuous-time systems. 

Obviously, (l32l) has the same pattern as (fill) , so 

P 1 {to\t-i)(P(to\t-i) — P(to\tN))P *(to\t-i) 

= i-(i+pHtoit-i^n-'gphtoit-!))- 1 

= V{I-[I + ^ 2 )~ 1 )V T , (33) 

where VS 2 V T is the singular value decomposition of P(to\t-i)^G T lZ~ l GP{to\t-\)^. 

Then, following the same steps as in Section[3] we can obtain the efficiency of observa¬ 
tion configurations for continuous-time systems. 

B The sensitivity of observation networks for continuous-time 

systems 

Consider the generalized continuous-time linear system: 

5x(t) = M(t,to)5x(to), 

with the corresponding forecast perturbation of observations evolving from dx(to ) 

5y(t ) = H(t)8x(t). 

To be brief, let Wo = P _1 (fi7-i) and W(t) = (see Appendix A), and define 

the magnitude of the perturbation of the initial state and observations respectively by 

||fe(to)||p-i(i 0 | t _ l) = (5x(t 0 ),-P _1 (to|t-i)5®(to)), 

pN 

INIjit-Rt)} = / 

■J to 
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Thus, the perturbation growth for continuous-time system can be measured by 





(H(t)5x(t), R 1 (t)H(t)Sx(t))dt 

(5x(to),fl 0 N M(t,t 0 ) T H(t) T R~ 1 (t)H(t)M(t,t 0 )6x(t 0 )dt) 
(6x(t 0 ),P- 1 (t 0 \t-i)Sx(t 0 )) 
{8x{t Q ),g*Ti- l g5x(tQ)) ^ 

(Sx(t 0 ),P-\t 0 \t. 1 )Sx(t 0 )) J ’ 


where Q and TZ~ 1 are defined in Appendix A. 

To find the directions maximizing the ratio, we need to find the solutions of the singular 
value problem: 


Pz{to\t-i)G T K 1 gP 2 (t Q \t- 1 )v k = s 2 k v k , 
gp-'itoit-i^uk = s 2 k u k . 


(34) 

(35) 


where si ^ S2 ^ ^ s n ^ 0, {v k \ r j l =l and {u k }™ =1 arc orthogonal singular vectors. 

Compared (l34l ) with (1331) in Appendix A, similar analysis and conclusions as Section [5] 
can be extended to continuous-time systems. 
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